perm filename SUBSLM.F4[SAT,LCS] blob sn#496778 filedate 1981-07-22 generic text, type T, neo UTF8
	SUBROUTINE SS
	COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
	COMMON/SLM/A(2,200),B(2,200),C(2,200),I(200),I1(200),
	1 I2(200),KT,TT
	DATA L1/1/,L2/2/,L0/0/
C	FIND MAX AND MIN POINTS AND SET SLOPES
	ICOUNT=0
710	DO 100 K=2,N-1
	KP=K-1
	KQ=K+1
	IF(Y(K).GT.Y(KP).AND.Y(K).GT.Y(KQ)) GO TO 20
	IF(Y(K).LT.Y(KP).AND.Y(K).LT.Y(KQ)) GO TO 20
	IF(X(K).GT.X(KP).AND.X(K).GT.X(KQ)) GO TO 30
	IF(X(K).LT.X(KP).AND.X(K).LT.X(KQ)) GO TO 30
	GO TO 100
20	S(K)=0
	ICOUNT=ICOUNT+1
	I(ICOUNT)=K
	GO TO 100
30	S(K)=1001
	ICOUNT=ICOUNT+1
	I(ICOUNT)=K
100	CONTINUE
	ICOUNT=ICOUNT+1
	I(ICOUNT)=N
C	FIND POINTS SUCH THAT X(K)=X(K+1) OR Y(K)=Y(K+1)
	DO 900 K=2,N-1
	IF(Y(K).EQ.Y(K+1)) I2(K)=1
	IF(X(K).EQ.X(K+1)) I2(K)=1
900	CONTINUE
C	SET THE SLOPES BETWEEN
	IA=1
	IX=0
110	IF(IX.EQ.ICOUNT) GO TO 200
	IX=IX+1
	IA=IA+1
120	IF(IA.EQ.I(IX)) GO TO 110
	KP=IA+1
	S(IA)=(Y(KP)-Y(IA-1))/(X(KP)-X(IA-1))
	IA=KP
	GO TO 120
200	CONTINUE
C	SET BEGIN AND END SLOPES
CC	SZ=1002001.
	K=1
	CALL S101
	CALL S102(X(1),Y(1),X(2),Y(2),S(1),S(2))
	IF(S(K).EQ.-1001.) S(K)=1001.
	IF(S(K+1).EQ.-1001.) S(K+1)=1001.
	K=N-1
	CALL S101
	CALL S102(X(N),Y(N),X(N-1),Y(N-1),S(N),S(N-1))
	IF(S(K).EQ.-1001.) S(K)=1001.
	IF(S(K+1).EQ.-1001.) S(K+1)=1001.
C	RESET THE SLOPES FOR STRAIGHT LINES
	DO 610 K=1,N-2
	U1=SLM2(L0)
	U2=SLM2(L1)
	IF(ABS(U1-U2).GT..0001) GO TO 610
	I1(K)=1
	I1(K+1)=1
	S(K)=U1
	S(K+1)=U1
	S(K+2)=U1
610	CONTINUE
C	ADD POINTS
	K1=N-1
	DO 300 K=1,K1
	IF(I1(K).EQ.1) GO TO 300
	IF(I2(K).EQ.1) GO TO 300
	FLG=0.
	IF(S(K).EQ.0..AND.S(K+1).EQ.0.) FLG=1.
	IF(S(K).EQ.1001..AND.S(K+1).EQ.1001.) FLG=2.
	IF(FLG.EQ.1..OR.FLG.EQ.2.) GO TO 205
CC	SZ=1002001.
	CALL S101
	KP=K+1
	U=(Y(K)-Y(KP)-S(K)*(X(K)-X(KP)))/(S(KP)-S(K))
	AX=X(KP)+U
	AY=Y(KP)+S(KP)*U
	XA=X(K)
1610	XB=X(KP)
	IF(XA.LE.XB) GO TO 202
	XA=X(KP)
	XB=X(K)
202	YA=Y(K)
	YB=Y(KP)
	IF(YA.LE.YB) GO TO 204
	YA=Y(KP)
	YB=Y(K)
204	CONTINUE
	IF(AX.GE.XA.AND.AX.LE.XB.AND.
     1	AY.GE.YA.AND.AY.LE.YB) GO TO 290
205	K1=K1+1
	DO 210 K2=K+1,N
	K3=N+K+1-K2
	KP=K3+1
	X(KP)=X(K3)
	Y(KP)=Y(K3)
	I1(KP)=I1(K3)
	I2(KP)=I2(K3)
210	S(KP)=S(K3)
	IF(FLG.EQ.1..OR.FLG.EQ.2.) GO TO 280
	N=N+1
	KP=K+1
	KQ=K+2
	XK=X(K)
	XQ=X(KQ)
	XZ=XK-XQ
C	Z0=X(K)**2-X(KQ)**2
	Z0=XK*XK+XQ*XQ
C	Z1=Y(K)-Y(KQ)-S(K)*(X(K)-X(KQ))-((X(K)+X(KQ))/2.
C     1	-X(K))*(S(K)-S(KQ))
	Z1=Y(K)-Y(KQ)-S(K)*XZ-(XK+XQ)/2.-XK*(S(K)-S(KQ))
C     	Z2=X(K)**3-X(KQ)**3-3*(X(K)-X(KQ))*X(K)**2
C     1	-1.5*(X(K)+X(KQ))*Z0+3*X(K)*Z0
     	Z2=XK**3-XQ**3-3*(XK-XQ)*XK*XK-1.5*(XK+XQ)*Z0+3*XK*Z0
	AZ=Z1/Z2
C	BZ=(S(K)-S(KQ)-3*(X(K)**2-X(KQ)**2)*AZ)/
C     1	(2*(X(K)-X(KQ)))
	BZ=(S(K)-S(KQ)-3*(XK*XK-XQ*XQ)*AZ)/(2*XZ)
C	CZ=S(K)-3*AZ*X(K)**2-2*BZ*X(K)
	CZ=S(K)-3*AZ*XK*XK-2*BZ*XK
C	DZ=Y(K)-AZ*X(K)**3-BZ*X(K)**2-CZ*X(K)
	DZ=Y(K)-AZ*XK**3-BZ*XK*XK-CZ*XK
	X(KP)=-BZ/(3.*AZ)
	Y(KP)=AZ*X(KP)**3+BZ*X(KP)**2+CZ*X(KP)+DZ
	S(KP)=3*AZ*X(KP)**2+2*BZ*X(KP)+CZ
	IF(Y(KP).LT.YB.AND.Y(KP).GT.YA) GO TO 290
	X(KP)=(XQ+XK)/2.
	Y(KP)=(Y(KQ)+Y(K))/2.
	S(KP)=0.
	GO TO 290
280	N=N+1
	IF(FLG.EQ.2.) GO TO 282
	S(K+1)=2.*(Y(K+2)-Y(K))/(X(K+2)-X(K))
	GO TO 284
282	S(K+1)=(Y(K+2)-Y(K))/(2.*(X(K+2)-X(K)))
284	X(K+1)=(X(K+2)+X(K))/2.
	Y(K+1)=(Y(K+2)+Y(K))/2.
290	IF(S(K).EQ.-1001.) S(K)=1001.
	IF(S(K+1).EQ.-1001.) S(K+1)=1001.
300	CONTINUE
C	FIT THE POINTS WITH N-1 CURVES
305	DO 400 K=1,N-1
	KP=K+1
	IF(I1(K).EQ.0) GO TO 310
	A(1,K)=0.
	A(2,K)=0.
	B(1,K)=X(KP)-X(K)
	B(2,K)=Y(KP)-Y(K)
	C(1,K)=X(K)
	C(2,K)=Y(K)
	GO TO 400
310	CALL S101
	B(1,K)=2*(Y(KP)-Y(K)-S(KP)*(X(KP)-X(K)))/
     1	(S(K)-S(KP))
	B(2,K)=S(K)*B(1,K)
	A(1,K)=X(KP)-X(K)-B(1,K)
	A(2,K)=Y(KP)-Y(K)-S(K)*B(1,K)
	C(1,K)=X(K)
	C(2,K)=Y(K)
400	CONTINUE
C	CALCULATE 512 POINTS
	M=N-1
	R=M/511.
	T=-R
	DO 500 L=1,511
	T=T+R
	KT=T+1
	KU=T
	TT=T-KU
	X1(L)=SLM3(L1)
500	Y1(L)=SLM3(L2)
	X1(512)=X(N)
	Y1(512)=Y(N)
600	RETURN
	END


C	GIVEN TWO POINTS (XX1,YY1), (XX2,YY2) AND A SLOPE SS2
C	THIS ROUTINE FINDS THE SLOPE SS1
	SUBROUTINE S102(XX1,YY1,XX2,YY2,SS1,SS2)
	LF=0
	IF(SS2.NE.0.) GO TO 5
	IF((XX1.LT.XX2.AND.YY1.LT.YY2).OR.
     1	(XX2.LT.XX1.AND.YY2.LT.YY1)) SS1=1001.
	IF((XX1.LT.XX2.AND.YY1.GT.YY2).OR.
     1	(XX2.LT.XX1.AND.YY2.GT.YY1)) SS1=-1001.
	GO TO 100
5	T=(YY1-YY2)/SS2
	X=XX2+T
	Y=YY1
	XMAX=XX2
	XMIN=XX1
	IF(XX2.GE.XX1) GO TO 10
	XMAX=XX1
	XMIN=XX2
10	IF(X.GT.XMIN.AND.X.LT.XMAX) GO TO 20
	T=XX1-XX2
	Y=YY2+SS2*T
	X=XX1
	LF=1
20	IF(LF.EQ.1) GO TO 30
	D2=XX1-X
	D1=X-XX2
	R=D2/D1
	XX=XX2
	YY=(YY2+R*YY1)/(1.+R)
	GO TO 40
30	D2=YY1-Y
	D1=Y-YY2
	R=D2/D1
	YY=YY2
	XX=(XX2+R*XX1)/(1.+R)
40	SS1=(YY-YY1)/(XX-XX1)
100	RETURN
	END

	SUBROUTINE S101
	COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
	DATA SZ/1002001./
	SS=S(K)**2
	SSS=S(K+1)**2
	IF(SS.NE.SZ.AND.SSS.NE.SZ)RETURN
	IF(SS.EQ.SZ.AND.X(K).LT.X(K+1).AND.
     1	Y(K).GT.Y(K+1)) S(K)=-1001.
	IF(S(K)**2.EQ.SZ.AND.X(K).GT.X(K+1).AND.
     1	Y(K+1).GT.Y(K)) S(K)=-1001.
	IF(SSS.EQ.SZ.AND.X(K).LT.X(K+1).AND.
     1	Y(K).GT.Y(K+1)) S(K+1)=-1001.
	IF(S(K+1)**2.EQ.SZ.AND.X(K).GT.X(K+1).AND.
     1	Y(K).LT.Y(K+1)) S(K+1)=-1001.
100	RETURN
	END

	FUNCTION SLM2(L)
	COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
	M=K+L
	J=M+1
	SLM2=(Y(J)-Y(M))/(X(J)-X(M))
	END

	FUNCTION SLM3(L)
	COMMON/SLM/A(2,200),B(2,200),C(2,200),I(200),I1(200),
	1 I2(200),KT,TT
	SLM3=A(L,KT)*TT**2+B(L,KT)*TT+C(L,KT)
	END